######################################################################################################
## This program plots the descriptive statistics of covariates ##
## Model refers to "Distribution Regression with Selection"
## Authors: Chernozhukov, V., I. Fernandez-Val and S. Luo
## Last Update: 10/13/2017

## Mean by sex-couple-year -- time trend



rm(list=ls());

#######################################
# Import data from a proper directory #
#######################################
mydata <- read.table("drs_fesdata.csv", header=TRUE, sep= ",");

attach(mydata);
mean.gender.married <- aggregate(cbind(work, lw, uwage, tubeninc0, beninc0, age, ageced, yrseduc), by=list(sex,couple,year), FUN=mean, na.rm=T);
mean.gender.work    <- aggregate(cbind(couple, tubeninc0, beninc0, age, ageced, yrseduc, numch, numchyt6), by=list(sex,work,year), FUN=mean, na.rm=T);
mean.gender         <- aggregate(cbind(work, lw, uwage, tubeninc0, beninc0, couple, age, ageced, yrseduc, numch, numchyt6), by=list(sex,year), FUN=mean, na.rm=T);

kid <- ifelse(numch==0,0,1);
mean.gender.married.kid <- aggregate(cbind(work, lw, uwage, tubeninc0, beninc0, age, ageced, yrseduc), 
                                     by=list(sex,couple,kid,year),FUN=mean, na.rm=T);



q90.gender <- aggregate(lw, by=list(sex,year), FUN=quantile, probs=0.9, na.rm=T);
q10.gender <- aggregate(lw, by=list(sex,year), FUN=quantile, probs=0.1, na.rm=T);


bornyear = sort(unique(cohort));
pct.14 = sapply(bornyear,function(i){sum(ageced==14&cohort==i)/sum(cohort==i)});
pct.15 = sapply(bornyear,function(i){sum(ageced==15&cohort==i)/sum(cohort==i)});
pct.16 = sapply(bornyear,function(i){sum(ageced==16&cohort==i)/sum(cohort==i)});

cohort.educ = cohort + ageced;
leaveyear = sort(unique(cohort.educ));
pct.14.educ = sapply(leaveyear,function(i){sum(ageced==14&cohort.educ==i)/sum(cohort.educ==i)});
pct.15.educ = sapply(leaveyear,function(i){sum(ageced==15&cohort.educ==i)/sum(cohort.educ==i)});
pct.16.educ = sapply(leaveyear,function(i){sum(ageced==16&cohort.educ==i)/sum(cohort.educ==i)});

detach(mydata);




pdf("Inequality.pdf");
matplot(c(1980:2010),
        cbind(mean.gender$lw[mean.gender$Group.1 == 1],
              mean.gender$lw[mean.gender$Group.1 == 2])[3:33,],
        col=c("blue", "red"), type="l",lty=1, xlab="Year", ylab="log(Hourly Wage)", main="Log of Usual Hourly Wage");
legend("bottomright", legend=c("Male","Female"), col=c("blue","red"),lty=1);

matplot(c(1980:2010),
        cbind(q90.gender$x[mean.gender$Group.1 == 1] - q10.gender$x[mean.gender$Group.1 == 1],
              q90.gender$x[mean.gender$Group.1 == 2] - q10.gender$x[mean.gender$Group.1 == 2])[3:33,],
        col=c("blue", "red"), type="l",lty=1, xlab="Year", ylab="log of 90/10 Percentile Ratio", main="Inequality within Gender Distribution");
legend("bottomright", legend=c("Male","Female"), col=c("blue","red"),lty=1);

matplot(c(1980:2010),
        cbind(mean.gender$work[mean.gender$Group.1 == 1],
              mean.gender$work[mean.gender$Group.1 == 2])[3:33,],
        col=c("blue", "red"), type="l",lty=1, xlab="Year", ylab="Employment Rate", main="Employment Rate");
legend("bottomright", legend=c("Male","Female"), col=c("blue","red"),lty=1);


dev.off();




pdf("Summary_Stats_Trend.pdf");

# Employment Rate (work) #
matplot(c(1978:2013),
        cbind(mean.gender$work[mean.gender$Group.1 == 1],
              mean.gender$work[mean.gender$Group.1 == 2]),
        col=c("blue", "red"), type="l",lty=1, xlab="Year", ylab="Employment Rate", main="Employment Rate");
legend("bottomright", legend=c("Male","Female"), col=c("blue","red"),lty=1);

matplot(c(1978:2013),
        cbind(mean.gender.married$work[mean.gender.married$Group.1 == 1 & mean.gender.married$Group.2 == 1],
              mean.gender.married$work[mean.gender.married$Group.1 == 1 & mean.gender.married$Group.2 == 0],
              mean.gender.married$work[mean.gender.married$Group.1 == 2 & mean.gender.married$Group.2 == 1],
              mean.gender.married$work[mean.gender.married$Group.1 == 2 & mean.gender.married$Group.2 == 0]),
        col=c("blue", "blue", "red", "red"), type="l", lty=c(1,2,1,2), xlab="Year", ylab="Employment Rate", 
        main="Employment Rate -- Married/Single");
legend("bottomright", legend=c("Married Male","Single Male","Married Female","Single Female"), 
       col=c("blue","blue","red","red"), lty=c(1,2,1,2));

matplot(c(1978:2013),
        cbind(mean.gender.married.kid$work[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$work[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 0],
              mean.gender.married.kid$work[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$work[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 0],
              mean.gender.married.kid$work[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$work[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 0],
              mean.gender.married.kid$work[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$work[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 0]),
        col=c("blue", "blue", "green", "green", "red", "red", "orange", "orange"), type="l", lty=c(1,2,1,2,1,2,1,2), 
        xlab="Year", ylab="Employment Rate", main="Employment Rate -- Married x Kids");
legend("bottomright", legend=c("Married Male w/ Kids","Married Male w/o Kids","Single Male w/ Kids","Single Male w/o Kids",
                               "Married Female w/ Kids","Married Female w/o Kids","Single Female w/ Kids","Single Female w/o Kids"), 
       col=c("blue", "blue", "green", "green", "red", "red", "orange", "orange"), lty=c(1,2,1,2,1,2,1,2));



# Hourly Wage #
matplot(c(1978:2013),
        cbind(mean.gender$uwage[mean.gender$Group.1 == 1],
              mean.gender$uwage[mean.gender$Group.1 == 2]),
        col=c("blue", "red"), type="l",lty=1, xlab="Year", ylab="Hourly Wage", main="Usual Hourly Wage");
legend("bottomright", legend=c("Male","Female"), col=c("blue","red"),lty=1);

matplot(c(1978:2013),
        cbind(mean.gender.married$uwage[mean.gender.married$Group.1 == 1 & mean.gender.married$Group.2 == 1],
              mean.gender.married$uwage[mean.gender.married$Group.1 == 1 & mean.gender.married$Group.2 == 0],
              mean.gender.married$uwage[mean.gender.married$Group.1 == 2 & mean.gender.married$Group.2 == 1],
              mean.gender.married$uwage[mean.gender.married$Group.1 == 2 & mean.gender.married$Group.2 == 0]),
        col=c("blue", "blue", "red", "red"), type="l", lty=c(1,2,1,2), xlab="Year", ylab="Hourly Wage", 
        main="Usual Hourly Wage -- Married/Single");
legend("bottomright", legend=c("Married Male","Single Male","Married Female","Single Female"), 
       col=c("blue","blue","red","red"), lty=c(1,2,1,2));


# Log of Hourly Wage #
matplot(c(1978:2013),
        cbind(mean.gender$lw[mean.gender$Group.1 == 1],
              mean.gender$lw[mean.gender$Group.1 == 2]),
        col=c("blue", "red"), type="l",lty=1, xlab="Year", ylab="log(Hourly Wage)", main="Log of Usual Hourly Wage");
legend("bottomright", legend=c("Male","Female"), col=c("blue","red"),lty=1);

matplot(c(1978:2013),
        cbind(mean.gender.married$lw[mean.gender.married$Group.1 == 1 & mean.gender.married$Group.2 == 1],
              mean.gender.married$lw[mean.gender.married$Group.1 == 1 & mean.gender.married$Group.2 == 0],
              mean.gender.married$lw[mean.gender.married$Group.1 == 2 & mean.gender.married$Group.2 == 1],
              mean.gender.married$lw[mean.gender.married$Group.1 == 2 & mean.gender.married$Group.2 == 0]),
        col=c("blue", "blue", "red", "red"), type="l", lty=c(1,2,1,2), xlab="Year", ylab="log(Hourly Wage)", 
        main="Log of Usual Hourly Wage -- Married/Single");
legend("bottomright", legend=c("Married Male","Single Male","Married Female","Single Female"), 
       col=c("blue","blue","red","red"), lty=c(1,2,1,2));

matplot(c(1978:2013),
        cbind(mean.gender.married.kid$lw[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$lw[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 0],
              mean.gender.married.kid$lw[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$lw[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 0],
              mean.gender.married.kid$lw[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$lw[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 0],
              mean.gender.married.kid$lw[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$lw[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 0]),
        col=c("blue", "blue", "green", "green", "red", "red", "orange", "orange"), type="l", lty=c(1,2,1,2,1,2,1,2), 
        xlab="Year", ylab="log(Hourly Wage)", main="Log of Usual Hourly Wage -- Married x Kids");
legend("bottomright", legend=c("Married Male w/ Kids","Married Male w/o Kids","Single Male w/ Kids","Single Male w/o Kids",
                               "Married Female w/ Kids","Married Female w/o Kids","Single Female w/ Kids","Single Female w/o Kids"), 
       col=c("blue", "blue", "green", "green", "red", "red", "orange", "orange"), lty=c(1,2,1,2,1,2,1,2));


# Marriage Rate #
matplot(c(1978:2013),
        cbind(mean.gender$couple[mean.gender$Group.1 == 1],
              mean.gender$couple[mean.gender$Group.1 == 2]),
        col=c("blue", "red"), type="l",lty=1, xlab="Year", ylab="Marriage Rate", main="Marriage Rate");
legend("bottomleft", legend=c("Male","Female"), col=c("blue","red"),lty=1);

matplot(c(1978:2013),
        cbind(mean.gender.work$couple[mean.gender.work$Group.1 == 1 & mean.gender.work$Group.2 == 1],
              mean.gender.work$couple[mean.gender.work$Group.1 == 1 & mean.gender.work$Group.2 == 0],
              mean.gender.work$couple[mean.gender.work$Group.1 == 2 & mean.gender.work$Group.2 == 1],
              mean.gender.work$couple[mean.gender.work$Group.1 == 2 & mean.gender.work$Group.2 == 0]),
        col=c("blue", "blue", "red", "red"), type="l", lty=c(1,2,1,2), xlab="Year", ylab="Marriage Rate", 
        main="Marriage Rate -- Employed/Nonemployed");
legend("bottomleft", legend=c("Employed Male","Nonemployed Male","Employed Female","Nonemployed Female"), 
       col=c("blue","blue","red","red"), lty=c(1,2,1,2));


# Number of Children #
matplot(c(1978:2013),
        cbind(mean.gender$numch[mean.gender$Group.1 == 1],
              mean.gender$numch[mean.gender$Group.1 == 2]),
        col=c("blue", "red"), type="l",lty=1, xlab="Year", ylab="# of Children", main="Average Number of Children");
legend("bottomleft", legend=c("Male","Female"), col=c("blue","red"),lty=1);

matplot(c(1978:2013),
        cbind(mean.gender.work$numch[mean.gender.work$Group.1 == 1 & mean.gender.work$Group.2 == 1],
              mean.gender.work$numch[mean.gender.work$Group.1 == 1 & mean.gender.work$Group.2 == 0],
              mean.gender.work$numch[mean.gender.work$Group.1 == 2 & mean.gender.work$Group.2 == 1],
              mean.gender.work$numch[mean.gender.work$Group.1 == 2 & mean.gender.work$Group.2 == 0]),
        col=c("blue", "blue", "red", "red"), type="l", lty=c(1,2,1,2), xlab="Year", ylab="# of Children", 
        main="Average Number of Children -- Employed/Nonemployed");
legend("bottomleft", legend=c("Employed Male","Nonemployed Male","Employed Female","Nonemployed Female"), 
       col=c("blue","blue","red","red"), lty=c(1,2,1,2));


# Number of Young Children #
matplot(c(1978:2013),
        cbind(mean.gender$numchyt6[mean.gender$Group.1 == 1],
              mean.gender$numchyt6[mean.gender$Group.1 == 2]),
        col=c("blue", "red"), type="l",lty=1, xlab="Year", ylab="# of Young Children", main="Number of Young Children(< 7)");
legend("bottomleft", legend=c("Male","Female"), col=c("blue","red"),lty=1);

matplot(c(1978:2013),
        cbind(mean.gender.work$numchyt6[mean.gender.work$Group.1 == 1 & mean.gender.work$Group.2 == 1],
              mean.gender.work$numchyt6[mean.gender.work$Group.1 == 1 & mean.gender.work$Group.2 == 0],
              mean.gender.work$numchyt6[mean.gender.work$Group.1 == 2 & mean.gender.work$Group.2 == 1],
              mean.gender.work$numchyt6[mean.gender.work$Group.1 == 2 & mean.gender.work$Group.2 == 0]),
        col=c("blue", "blue", "red", "red"), type="l", lty=c(1,2,1,2), xlab="Year", ylab="# of Young Children", 
        main="Number of Young Children(< 7) -- Employed/Nonemployed");
legend("bottomleft", legend=c("Employed Male","Nonemployed Male","Employed Female","Nonemployed Female"), 
       col=c("blue","blue","red","red"), lty=c(1,2,1,2));



# Out-of-Work Income #
matplot(c(1978:2013),
        cbind(mean.gender$tubeninc0[mean.gender$Group.1 == 1],
              mean.gender$tubeninc0[mean.gender$Group.1 == 2]),
        col=c("blue", "red"), type="l",lty=1, xlab="Year", ylab="Benefit Income", 
        main="Out-of-Work Income within Tax Unit");
legend("bottomright", legend=c("Male","Female"), col=c("blue","red"),lty=1);

matplot(c(1978:2013),
        cbind(mean.gender.married$tubeninc0[mean.gender.married$Group.1 == 1 & mean.gender.married$Group.2 == 1],
              mean.gender.married$tubeninc0[mean.gender.married$Group.1 == 1 & mean.gender.married$Group.2 == 0],
              mean.gender.married$tubeninc0[mean.gender.married$Group.1 == 2 & mean.gender.married$Group.2 == 1],
              mean.gender.married$tubeninc0[mean.gender.married$Group.1 == 2 & mean.gender.married$Group.2 == 0]),
        col=c("blue", "blue", "red", "red"), type="l", lty=c(1,2,1,2), xlab="Year", ylab="Benefit Income", 
        main="Out-of-Work Income within Tax Unit -- Married/Single");
legend("bottomright", legend=c("Married Male","Single Male","Married Female","Single Female"), 
       col=c("blue","blue","red","red"), lty=c(1,2,1,2));

matplot(c(1978:2013),
        cbind(mean.gender.married.kid$tubeninc0[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$tubeninc0[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 0],
              mean.gender.married.kid$tubeninc0[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$tubeninc0[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 0],
              mean.gender.married.kid$tubeninc0[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$tubeninc0[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 0],
              mean.gender.married.kid$tubeninc0[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$tubeninc0[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 0]),
        col=c("blue", "blue", "green", "green", "red", "red", "orange", "orange"), type="l", lty=c(1,2,1,2,1,2,1,2), 
        xlab="Year", ylab="Benefit Income", main="Out-of-Work Income within Tax Unit -- Married x Kids");
legend("bottomright", legend=c("Married Male w/ Kids","Married Male w/o Kids","Single Male w/ Kids","Single Male w/o Kids",
                               "Married Female w/ Kids","Married Female w/o Kids","Single Female w/ Kids","Single Female w/o Kids"), 
       col=c("blue", "blue", "green", "green", "red", "red", "orange", "orange"), lty=c(1,2,1,2,1,2,1,2));


matplot(c(1978:2013),
        cbind(mean.gender.work$tubeninc0[mean.gender.work$Group.1 == 1 & mean.gender.work$Group.2 == 1],
              mean.gender.work$tubeninc0[mean.gender.work$Group.1 == 1 & mean.gender.work$Group.2 == 0],
              mean.gender.work$tubeninc0[mean.gender.work$Group.1 == 2 & mean.gender.work$Group.2 == 1],
              mean.gender.work$tubeninc0[mean.gender.work$Group.1 == 2 & mean.gender.work$Group.2 == 0]),
        col=c("blue", "blue", "red", "red"), type="l", lty=c(1,2,1,2), xlab="Year", ylab="Benefit Income", 
        main="Out-of-Work Income within Tax Unit -- Employed/Nonemployed");
legend("bottomright", legend=c("Employed Male","Nonemployed Male","Employed Female","Nonemployed Female"), 
       col=c("blue","blue","red","red"), lty=c(1,2,1,2));



# Individual Out-of-Work Income #
matplot(c(1978:2013),
        cbind(mean.gender$beninc0[mean.gender$Group.1 == 1],
              mean.gender$beninc0[mean.gender$Group.1 == 2]),
        col=c("blue", "red"), type="l",lty=1, xlab="Year", ylab="Benefit Income", 
        main="Individual Out-of-Work Income");
legend("bottomright", legend=c("Male","Female"), col=c("blue","red"),lty=1);

matplot(c(1978:2013),
        cbind(mean.gender.married$beninc0[mean.gender.married$Group.1 == 1 & mean.gender.married$Group.2 == 1],
              mean.gender.married$beninc0[mean.gender.married$Group.1 == 1 & mean.gender.married$Group.2 == 0],
              mean.gender.married$beninc0[mean.gender.married$Group.1 == 2 & mean.gender.married$Group.2 == 1],
              mean.gender.married$beninc0[mean.gender.married$Group.1 == 2 & mean.gender.married$Group.2 == 0]),
        col=c("blue", "blue", "red", "red"), type="l", lty=c(1,2,1,2), xlab="Year", ylab="Benefit Income", 
        main="Individual Out-of-Work Income -- Married/Single");
legend("bottomright", legend=c("Married Male","Single Male","Married Female","Single Female"), 
       col=c("blue","blue","red","red"), lty=c(1,2,1,2));

matplot(c(1978:2013),
        cbind(mean.gender.married.kid$beninc0[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$beninc0[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 0],
              mean.gender.married.kid$beninc0[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$beninc0[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 0],
              mean.gender.married.kid$beninc0[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$beninc0[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 0],
              mean.gender.married.kid$beninc0[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$beninc0[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 0]),
        col=c("blue", "blue", "green", "green", "red", "red", "orange", "orange"), type="l", lty=c(1,2,1,2,1,2,1,2), 
        xlab="Year", ylab="Benefit Income", main="Individual Out-of-Work Income -- Married x Kids");
legend("bottomright", legend=c("Married Male w/ Kids","Married Male w/o Kids","Single Male w/ Kids","Single Male w/o Kids",
                               "Married Female w/ Kids","Married Female w/o Kids","Single Female w/ Kids","Single Female w/o Kids"), 
       col=c("blue", "blue", "green", "green", "red", "red", "orange", "orange"), lty=c(1,2,1,2,1,2,1,2));

matplot(c(1978:2013),
        cbind(mean.gender.work$beninc0[mean.gender.work$Group.1 == 1 & mean.gender.work$Group.2 == 1],
              mean.gender.work$beninc0[mean.gender.work$Group.1 == 1 & mean.gender.work$Group.2 == 0],
              mean.gender.work$beninc0[mean.gender.work$Group.1 == 2 & mean.gender.work$Group.2 == 1],
              mean.gender.work$beninc0[mean.gender.work$Group.1 == 2 & mean.gender.work$Group.2 == 0]),
        col=c("blue", "blue", "red", "red"), type="l", lty=c(1,2,1,2), xlab="Year", ylab="Benefit Income", 
        main="Individual Out-of-Work Income -- Employed/Nonemployed");
legend("bottomright", legend=c("Employed Male","Nonemployed Male","Employed Female","Nonemployed Female"), 
       col=c("blue","blue","red","red"), lty=c(1,2,1,2));




# Age #
matplot(c(1978:2013),
        cbind(mean.gender$age[mean.gender$Group.1 == 1],
              mean.gender$age[mean.gender$Group.1 == 2]),
        col=c("blue", "red"), type="l",lty=1, xlab="Year", ylab="Age", 
        main="Average Age");
legend("bottomright", legend=c("Male","Female"), col=c("blue","red"),lty=1);

matplot(c(1978:2013),
        cbind(mean.gender.married$age[mean.gender.married$Group.1 == 1 & mean.gender.married$Group.2 == 1],
              mean.gender.married$age[mean.gender.married$Group.1 == 1 & mean.gender.married$Group.2 == 0],
              mean.gender.married$age[mean.gender.married$Group.1 == 2 & mean.gender.married$Group.2 == 1],
              mean.gender.married$age[mean.gender.married$Group.1 == 2 & mean.gender.married$Group.2 == 0]),
        col=c("blue", "blue", "red", "red"), type="l", lty=c(1,2,1,2), xlab="Year", ylab="Age", 
        main="Average Age -- Married/Single");
legend("bottomright", legend=c("Married Male","Single Male","Married Female","Single Female"), 
       col=c("blue","blue","red","red"), lty=c(1,2,1,2));

matplot(c(1978:2013),
        cbind(mean.gender.married.kid$age[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$age[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 0],
              mean.gender.married.kid$age[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$age[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 0],
              mean.gender.married.kid$age[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$age[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 0],
              mean.gender.married.kid$age[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$age[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 0]),
        col=c("blue", "blue", "green", "green", "red", "red", "orange", "orange"), type="l", lty=c(1,2,1,2,1,2,1,2), 
        xlab="Year", ylab="Age", main="Average Age -- Married x Kids");
legend("bottomright", legend=c("Married Male w/ Kids","Married Male w/o Kids","Single Male w/ Kids","Single Male w/o Kids",
                               "Married Female w/ Kids","Married Female w/o Kids","Single Female w/ Kids","Single Female w/o Kids"), 
       col=c("blue", "blue", "green", "green", "red", "red", "orange", "orange"), lty=c(1,2,1,2,1,2,1,2));

matplot(c(1978:2013),
        cbind(mean.gender.work$age[mean.gender.work$Group.1 == 1 & mean.gender.work$Group.2 == 1],
              mean.gender.work$age[mean.gender.work$Group.1 == 1 & mean.gender.work$Group.2 == 0],
              mean.gender.work$age[mean.gender.work$Group.1 == 2 & mean.gender.work$Group.2 == 1],
              mean.gender.work$age[mean.gender.work$Group.1 == 2 & mean.gender.work$Group.2 == 0]),
        col=c("blue", "blue", "red", "red"), type="l", lty=c(1,2,1,2), xlab="Year", ylab="Age", 
        main="Average Age -- Employed/Nonemployed");
legend("bottomright", legend=c("Employed Male","Nonemployed Male","Employed Female","Nonemployed Female"), 
       col=c("blue","blue","red","red"), lty=c(1,2,1,2));


# Education #
matplot(c(1978:2013),
        cbind(mean.gender$ageced[mean.gender$Group.1 == 1],
              mean.gender$ageced[mean.gender$Group.1 == 2]),
        col=c("blue", "red"), type="l",lty=1, xlab="Year", ylab="Age Ceased Education", 
        main="Age Ceased Education");
legend("bottomright", legend=c("Male","Female"), col=c("blue","red"),lty=1);

matplot(c(1978:2013),
        cbind(mean.gender.married$ageced[mean.gender.married$Group.1 == 1 & mean.gender.married$Group.2 == 1],
              mean.gender.married$ageced[mean.gender.married$Group.1 == 1 & mean.gender.married$Group.2 == 0],
              mean.gender.married$ageced[mean.gender.married$Group.1 == 2 & mean.gender.married$Group.2 == 1],
              mean.gender.married$ageced[mean.gender.married$Group.1 == 2 & mean.gender.married$Group.2 == 0]),
        col=c("blue", "blue", "red", "red"), type="l", lty=c(1,2,1,2), xlab="Year", ylab="Age Ceased Education", 
        main="Age Ceased Education -- Married/Single");
legend("bottomright", legend=c("Married Male","Single Male","Married Female","Single Female"), 
       col=c("blue","blue","red","red"), lty=c(1,2,1,2));

matplot(c(1978:2013),
        cbind(mean.gender.married.kid$ageced[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$ageced[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 0],
              mean.gender.married.kid$ageced[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$ageced[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 0],
              mean.gender.married.kid$ageced[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$ageced[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 0],
              mean.gender.married.kid$ageced[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$ageced[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 0]),
        col=c("blue", "blue", "green", "green", "red", "red", "orange", "orange"), type="l", lty=c(1,2,1,2,1,2,1,2), 
        xlab="Year", ylab="Age Ceased Education", main="Age Ceased Education -- Married x Kids");
legend("bottomright", legend=c("Married Male w/ Kids","Married Male w/o Kids","Single Male w/ Kids","Single Male w/o Kids",
                               "Married Female w/ Kids","Married Female w/o Kids","Single Female w/ Kids","Single Female w/o Kids"), 
       col=c("blue", "blue", "green", "green", "red", "red", "orange", "orange"), lty=c(1,2,1,2,1,2,1,2));

matplot(c(1978:2013),
        cbind(mean.gender.work$ageced[mean.gender.work$Group.1 == 1 & mean.gender.work$Group.2 == 1],
              mean.gender.work$ageced[mean.gender.work$Group.1 == 1 & mean.gender.work$Group.2 == 0],
              mean.gender.work$ageced[mean.gender.work$Group.1 == 2 & mean.gender.work$Group.2 == 1],
              mean.gender.work$ageced[mean.gender.work$Group.1 == 2 & mean.gender.work$Group.2 == 0]),
        col=c("blue", "blue", "red", "red"), type="l", lty=c(1,2,1,2), xlab="Year", ylab="Age Ceased Education", 
        main="Age Ceased Education -- Employed/Nonemployed");
legend("bottomright", legend=c("Employed Male","Nonemployed Male","Employed Female","Nonemployed Female"), 
       col=c("blue","blue","red","red"), lty=c(1,2,1,2));



# Years of Schooling #
matplot(c(1978:2013),
        cbind(mean.gender$yrseduc[mean.gender$Group.1 == 1],
              mean.gender$yrseduc[mean.gender$Group.1 == 2]),
        col=c("blue", "red"), type="l",lty=1, xlab="Year", ylab="Years of Schooling", 
        main="Years of Schooling");
legend("bottomright", legend=c("Male","Female"), col=c("blue","red"),lty=1);

matplot(c(1978:2013),
        cbind(mean.gender.married$yrseduc[mean.gender.married$Group.1 == 1 & mean.gender.married$Group.2 == 1],
              mean.gender.married$yrseduc[mean.gender.married$Group.1 == 1 & mean.gender.married$Group.2 == 0],
              mean.gender.married$yrseduc[mean.gender.married$Group.1 == 2 & mean.gender.married$Group.2 == 1],
              mean.gender.married$yrseduc[mean.gender.married$Group.1 == 2 & mean.gender.married$Group.2 == 0]),
        col=c("blue", "blue", "red", "red"), type="l", lty=c(1,2,1,2), xlab="Year", ylab="Years of Schooling", 
        main="Years of Schooling -- Married/Single");
legend("bottomright", legend=c("Married Male","Single Male","Married Female","Single Female"), 
       col=c("blue","blue","red","red"), lty=c(1,2,1,2));

matplot(c(1978:2013),
        cbind(mean.gender.married.kid$yrseduc[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$yrseduc[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 0],
              mean.gender.married.kid$yrseduc[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$yrseduc[mean.gender.married.kid$Group.1 == 1 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 0],
              mean.gender.married.kid$yrseduc[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$yrseduc[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 1 & mean.gender.married.kid$Group.3 == 0],
              mean.gender.married.kid$yrseduc[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 1],
              mean.gender.married.kid$yrseduc[mean.gender.married.kid$Group.1 == 2 & mean.gender.married.kid$Group.2 == 0 & mean.gender.married.kid$Group.3 == 0]),
        col=c("blue", "blue", "green", "green", "red", "red", "orange", "orange"), type="l", lty=c(1,2,1,2,1,2,1,2), 
        xlab="Year", ylab="Years of Schooling", main="Years of Schooling -- Married x Kids");
legend("bottomright", legend=c("Married Male w/ Kids","Married Male w/o Kids","Single Male w/ Kids","Single Male w/o Kids",
                               "Married Female w/ Kids","Married Female w/o Kids","Single Female w/ Kids","Single Female w/o Kids"), 
       col=c("blue", "blue", "green", "green", "red", "red", "orange", "orange"), lty=c(1,2,1,2,1,2,1,2));


matplot(c(1978:2013),
        cbind(mean.gender.work$yrseduc[mean.gender.work$Group.1 == 1 & mean.gender.work$Group.2 == 1],
              mean.gender.work$yrseduc[mean.gender.work$Group.1 == 1 & mean.gender.work$Group.2 == 0],
              mean.gender.work$yrseduc[mean.gender.work$Group.1 == 2 & mean.gender.work$Group.2 == 1],
              mean.gender.work$yrseduc[mean.gender.work$Group.1 == 2 & mean.gender.work$Group.2 == 0]),
        col=c("blue", "blue", "red", "red"), type="l", lty=c(1,2,1,2), xlab="Year", ylab="Years of Schooling", 
        main="Years of Schooling -- Employed/Nonemployed");
legend("bottomright", legend=c("Employed Male","Nonemployed Male","Employed Female","Nonemployed Female"), 
       col=c("blue","blue","red","red"), lty=c(1,2,1,2));



dev.off();



pdf("LeavingSchool.pdf");
  matplot(bornyear, cbind(pct.14,pct.15,pct.16)*100, type="l", lty=1,
          xlab="Cohort", ylab="% of Cohort", main="Percentage of Subpopulation within Birth Cohort");
  legend("topright",c("Leaving at 14","Leaving at 15","Leaving at 16"), col=c(1:3),lty=1);

  matplot(leaveyear, cbind(pct.14.educ,pct.15.educ,pct.16.educ)*100, type="l", lty=1, col=c(2:4), xlim=c(1930,2010), 
         xlab="Year of Leaving School", ylab="% of Cohort", main="Percentage of Subpopulation within Education Cohort");
  legend("topright",c("Leaving at 14","Leaving at 15","Leaving at 16"), col=c(2:4),lty=1);
  abline(v=1947,lty=2);
  abline(v=1972,lty=2);
  axis(side=1, at=c(1947,1972));
dev.off();
